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1. What is CosMoPMC? 



1 . What is CosMoPMC? 



inamei 

CosMoPMC (Cosmology Population Monte Carlo) is a Bayesian sampling method to explore 
the likelihood of various cosmological probes. The sampling engine is implemented with the 
package pmclib. It is called Population Monte Carlo (PMC), which is a novel technique to sample 
from the posterior ( Cappe et al.|2008 ). PMC is an adaptive importance sampling method which 
iteratively improves the proposal to approximate the posterior. This code has been introduced, 



tested and applied to various cosmology data sets in Wraith et al. (2009 1. Results on the Bayesian 



evidence using PMC are discussed in Kilbinger et al. (2010 1. 



1.1. Importance sampling 

One of the main goals in Bayesian inference is to obtain integrals of the form 

n(f) = J f(x)7:(x)dx 



(1) 



over the posterior distribution k which depends on the /^-dimensional parameter x, where / is 
an arbitrary function with finite expectation under n. Of interest are for example the parameter 
mean (/ = id) or confidence regions S with / = I5 being the indicator function of S . The 
Bayesian evidence E, used in model comparison techniques, is obtained by setting f = I, but 
instead of n using the unnormalised posterior n' = L • P in ([l]), with L being the likelihood and 
P the prior. 

The evaluation of ([T]l is challenging because the posterior is in general not available analytically, 
and the parameter space can be high-dimensional. Monte-Carlo methods to approximate the 
above integrals consist in providing a sample {xn}n=i...N under n, and approximating ([T]l by the 
estimator 



1 

n=l 



(2) 



Markov Chain Monte Carlo (MCMC) produces a Markov chain of points for which n is the 



limiting distribution. The popular and widely-used package cosmomc (http : //cosmologist . 
[info/cosmomc; Lewis & Bridle [2002[ ) implements MCMC exploration of the cosmological 
parameter space. 

Importance sampling on the other hand uses the identity 



J 



^(/) - /(x)7r(x)dx 



r 7t(x) 



q{x) 



(3) 
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where q is any probability density function with support including the support of n. A sample 
{xn] under q is then used to obtain the estimator 



N 

= (4) 

n=l 



1 v-i T^{x.n) 
^(/) ^ Tt / J(Xn)Wn; Wn ^ — — 



The function q is called the proposal or importance function, the quantities w„ are the importance 
weights. Population Monte Carlo (PMC) produces a sequence q' of importance functions (f = 



1 . . . r) to approximate the posterior n. Details of this algorithm are discussed in Wraith et al. 
( [20091 ). 



The package CosmoPMC provides a C-code for sampling and exploring the cosmological pa- 
rameter space using Population Monte Carlo. The code uses MPI to parallelize the calculation 
of the likelihood function. There is very little overhead and on a massive cluster the reduc- 
tion in wall-clock time can be enormous. Included in the package are post-processing, plotting 
and various other analysis scripts and programs. It also provides a Markov Chain Monte-Carlo 
sampler. 



1.2. This manual 

This manual describes the code CosmoPMC, and can be obtained from |www . cosmopmc . info 
CosmoPMC is the cosmology interface to the Population Monte Carlo (PMC) engine pmclib. 
Documentation on the PMC library can be found at the same url. The cosmology module of 
CosmoPMC can be used as stand-alone program, it has the name nicaea (http : //www2 . iap7| 
fr/users/kilbinge/ nicaea ). 

Warning: Use undocumented features of the code at your own risk! 



2. Installing CosmoPMC 
2.1. Software requirements 

CosmoPMC has been developed on GNU/Linux and Darwin/FreeBSD systems and should run 
on those architectures. Required are: 

• C-compiler (e.g. gcc, ice) 

• PMCLIB (Sect. |2.2| ) 

• GSL ( |http : //www . gnu . org/sof tware/gsT) , version 1.15 or higher 



FFTW dhttp : //www . f f tw . org |l 

Message Parsing Interface (MPI) ( jhttp : / /www-unix . itics . anl . gov/mpi ) for parallel 



calculations 
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2. Jnstaiiing CosMoPMC 



Optional: 

• csh, for post-processing, auxiliary scripts; recommended 



perl ( |http : //www . perl . orgp , for post-processing, auxiliary scripts; recommended 



yorick ( [http://yorick. sourceforge .net] ), post-processing, mainly plotting 



python (http://www.python.orgi, for running the configuration script 



R (http://www.r-project.org), post-processing 



To produce ID and 2D marginal posterior plots with scripts that come with CosmoPMC, either 
yorick or R are required. 

Necessary for CMB anisotropics support: 

• Fortran compiler (e.g. ifort) 

• Intel Math Kernel libraries ( |http : //software . Intel . com/en-us/intel-inkl| ) 



CAMB Chttp : //camb . info , http : //cosmologist . inf o/cosmoitic]) 



WMAP data and likehhood code (http : //lambda . gsf c . nasa . gov| ) 



2.2. Download and install pmclib 



The package pmclib can be downloaded from the CosmoPMC site http: //www. cosmopmcj 
Unfo 

After downloading, unpack the gzipped tar archive 

> tar xzf pmclib_x.y.tar.gz 

This creates the pmclib root directory pmclib_x . y. pmclib uses wai^instead of configure/make 
to compile and build the software. Change to that directory and type 

> ./waf --local configure 

See . /waf --help for options. The packages lua, hdf 5 and lapack are optionally linked with 
pmclib but are not necessary to run CosmoPMC. Corresponding warnings of missing files can be 
ignored. Instead of a local installation (indicated by --local), a install prefix can be specified 
with --prefix=PREFIX (default /usr/local). 



http : //code . google . com/p/waf 
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2. Jnstaiiing CosMoPMC 



2.3. Patch PMCLiB 

For CosMoPMC vl.2 and pmclib vl.x, a patch of the latter is necessary. From http : //www. | 
[cosmopmc . inf o| , download patch_pmclib_l . x_l . 2 . tar . gz and follow the instructions in 
the readme file readme_patch_pmclib_l . x_l . 2 . txt. 



2.4. Download and install CosmoPMC 

The newest version of CosmoPMC can be downloaded from the site |http : / /www . cosmopmc . 
linfo 

First, unpack the gzipped tar archive 

> tar xzf CosmoPMC.vl . 2 . tar . gz 

This creates the the CosmoPMC root directory CosmoPMC_vl . 2. Change to that directory and 
run 

> [python] ./configure.py 

This (poor man's) configure script copies the file Makefile .no Jiost to Makefile .host and 
sets host-specific variables and flags as given by the command-hne arguments. For a complete 
list, see 'configure.py --help'. 

Alternatively, you can copy by hand the file Makefile .no_host to Makefile .host and edit 
it. If the flags in this file are not sufficient to successfully compile the code, you can add more 
flags by rerunning configure . py, or by manually editing Makefile . main. Note that a flag in 
Makefile .main is overwritten if the same flag is present in Makefile .host. 

To compile the code, run 

> make; make clean 

On success, symbohc links to the binary executables (in . /exec) will be set in . /bin. 

It is convenient to define the environment variable COSMOPMC and to set it to the main CosmoPMC 
directory. For example, in the C-shell: 

> setenv COSMOPMC /path/to/CosmoPMC_vl . 2 

This command can be placed into the startup file (e.g. " / . cshrc for the C-shell). One can also 
add $COSMOPMC/bin to the PATH environment variable. 
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3. Running CosmoPMC 

3.1. Quick reference guide 
Examples 

To get familiar with CosmoPMC, use the examples which are contained in the package. Simply 
change to one of the subdirectories in $COSMOPMC/Demo/MCJDemo and proceed on to the point 
Run below. 



User-defined runs 

To run different likelihood combinations, or your own data, the following two steps are necessary 
to set up a CosmoPMC run. 

1. Data and parameter files 



Create new directory with newdir_pmc . sh When asked, enter the likelihood/data type. 



More than one type can be chosen by adding the corresponding (bit-coded) type id's. 
Symbolic links to corresponding files in $COSMOPMC/data are set, and parameter files 
from $COSMOPMC/par .files are copied to the new directory on request. 

If necessary, copy different or additional data and/or parameter files to the present direc- 
tory. 

2. Configuration file 

Create the PMC configuration file config_pmc. Examples for existing data modules can 
be found in $COSMOPMC/Demo/MC_Demo, see also Sect.[5]for details. 

In some cases, information about the galaxy redshift distribution(s) have to be provided, 
and the corresponding files copied (see $COSMOPMC/Demo for example files 'nofz*'). 



Run 

Type 

> $COSMOPMC/bin/cosmo _pmc.pl -n NCPU 



to run CosmoPMC on NCPU CPUs. See cosmo_pmc.pl -h' for more options. Depending on 



the type of initial proposal (Sect. 3.2 1, a maximum-search is started followed by a Fisher matrix 
calculation. After that, PMC is started. Fig.[T]shows a flow chart of the script's actions. 
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3. Running CosMoPMC 



Diagnostics 



Check the files perplexity and enc. If the perplexity reaches values of 0.8 or larger, and if the 
effective number of components (ENC) is not smaller than 1.5, the posterior has very likely been 
explored sufficiently. Those and other files are updated during run-time and can be monitored 
while PMC is running. See Sect. 3.3.1 for more details. 



Results 



The text file iter_|niter-l)/mean contains mean and confidence levels. The file 
iter_|niter-l}/all_contour2d.pdf shows the Id- and 2d-marginals. Plots can be redone 
or refined, or created from other than the last iteration with |plot_contour2d . pl[ Note that in 
the default setting, the posterior plots are not smoothed. See Sect. |6.1.T]For more details, and for 



information on the alternative script plot_conf idence . R 



3.2. CosMoPMC in detail 



This section describes in more detail how PMC is run, and which decisions the user has to make 
before starting and after stopping a PMC run. 



Initial proposal The choice of the initial proposal, used during the first PMC iteration, is of 
great importance for a successful PMC run. The following options are implemented, determined 
by the key 'sinitial' in the configuration file (see Sect.|5]l: 

1. sinitial = fisher jrshift The Fisher matrix is used as the covariance of a multi- 
variate Gaussian/Student-f distribution g. A mixture-model is constructed by creating D 
copies of g. Each copy is displaced from the ML point by a random uniform shift, and its 
variance is stretched by random uniform factor. 

2. sinitial = f isher_eigen A mixture-model is constructed in a similar way as the first 
case, with the difference that the shift from the ML point is now performed along the 
major axes of the Fisher ellipsoid. Note that if the Fisher matrix is diagonal, the shift of 
each component only concerns one parameter. 

3. sinitial = file The initial proposal is read from a file (of mixjnvdens format), e.g. 
from a previous PMC run. 

4. sinitial = random_pos Mixture-model components with random variance (up to half 
the box size) and random positions. This case should only be used if the posterior is 
suspected to be multi-modal, or the calculation of the Fisher matrix fails. 

In many cases, a mixture of multi-variate Gaussians as the proposal is the best choice. For that, 
set the degrees-of-freedom (y) parameter df to -1. For a posterior with heavy tails, a Student-f 
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Figure 1: Flow chart for cosmo_pmc .pi. 



distribution might be more suited. The degrees of freedom v can be chosen freely; y = 3 is a 
common choice. For v — > oo, a Gaussian distribution is reached asymptotically. 



If the Fisher matrix has to be calculated for the initial proposal, the script cosmo_pmc . pi calls 



max.post and go.fishing to estimate the maximum-likelihood point and the Hessian at that 



point, respectively. The script config_pmc_to_max_and_fish.pl can be used to create the 
corresponding configuration files from the PMC config file for manual calls of max_post and 
go_fishing. 



Updating the proposal The PMC algorithm automatically updates the proposal after each 
iteration, no user interference is necessary. 
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3. Running CosMoPMC 



The method to update the proposal is a variant of the Expectation-Maximization algorithm (EM, 
Dempster et al.|1977] l. It leads to an increase of the perplexity and an increase of ESS. Detailed 
descriptions of this algorithm in the case of multi-variate Gaussian and Student-f distributions 



can be found in Cappe et al. (2008) and Wraith et al. (20091. 



Dead components A component can 'die' during the updating if the number of points sam- 
pled from that component is less than MINCOUNT = 20, or its weight is smaller than the inverse 
total number of sample points 1 /N. There are two possibilities to proceed. First, the component 
is 'buried', its weight set to zero so that no points are sampled from it in subsequent iterations. 
Alternatively, the component can be revived. In this case, it is placed near the component (pdo 
which has maximum weight, and it is given the same covariance as (p^g . 

The first case is the standard method used in Wraith et al.| ( |2009] ). The second method tries to 
cure cases where the majority of components die. This can happen if they start too far off from 
the high-density posterior region. Often, only one component remains to the end, not capable of 
sampling the posterior reliably. 

Both options can be chosen using the config file (Sect.|5]l key sdead_comp = {bury | revive). 



Errors If an error occurs during the calculation of the likelihood, the error is intercepted and 
the likelihood is set to zero. Thus, the parameter vector for which the error occurs is attributed a 
zero importance weight and does not contribute to the final sample. An error message is printed 
to stderr (unless CosmoPMC is run with the option -q) and PMC continues with the next point. 

An error can be due to cosmological reasons, e.g. a redshift is probed which is larger than the 
maximum redshift in a loitering Universe. Further, a parameter could be outside the range of a 
fitting formulae, e.g. a very small scalar spectral index in the dark matter transfer function. 

Usually, the errors printed to stderr during PMC sampling can be ignored. 

Random numbers The GSL random number generator is used to generate random variables. 
It is initialised with a seed reading the current time, to produce different (pseudo-) random 
numbers at each call. The seed is written to the log file. Using the option '-s SEED', a user- 
specified seed can be defined. This is helpful if a run is to be repeated with identical results. 

3.3. Output files 

Each iteration / produces a number of output files which are stored in subdirectories iter_i of 
the CosmoPMC starting directory. Files which are not specific to a single iteration are placed in 
the starting directory. 
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3.3.1. Diagnostics 

Unlike in MCMC, with adaptive importance sampling one does not have to worry about con- 
vergence. In principle, the updating process can be stopped at any time. There are however 
diagnostics to indicate the quality and effectiveness of the sampling. 



Perplexity and effective sample size perplexity 



The perplexity p is defined in eq. (18) of Wraith et al. (20091. The range of p is [0; 1], and 



will approach unity if the proposal and posterior distribution are close together, as measured by 
the KuUback-Leibler divergence. The initial perplexity is typically very low (< 0.1) and should 
increase from iteration to iteration. Final values of 0.99 and larger are not uncommon, but also 
for p of about 0.6-0.8 very accurate results can be obtained. If p is smaller than say 0.1, the 
PMC sample is most likely not representative of the posterior. Intermediate values for p are not 
straight-forward to interpret. 

Closely related to the perplexity is the effective sample size ESS, which lies in the range [1; A^]. 
It is interpreted as the number of sample point with zero weight (Liu & Chen '1995). A large 
perplexity is usually accompanied by a high ESS. For a successful PMC run, ESS is much higher 
than the acceptance rate of a Monte Carlo Markov chain, which is typically between 0.15 and 
0.25. 

The file perplexity contains the iteration /, perplexity p, ESS for that iteration, and the total 
ESS. This file is updated after each iteration and can therefore be used to monitor a PMC run. 

If there are points with very large weights, they can dominate the other points whose normalised 
weights will be small. Even a few sample points might dominate the sum over weights and result 
in a low perplexity. The perplexity is the most sensitive quantity to those high-weight points, 
much more than e.g. the mean, the confidence intervals or the evidence. 



Effective number of proposal components enc 

The proposal q' provides useful information about the performance of a PMC run. For example, 
the effective number of components, defined in complete analogy to ESS, 



ENC 



d=i 



(5) 



is an indication of components with non-zero weight. If ENC is close to unity, the number of 
remaining components to sample the posterior is likely to be too small to provide a represen- 
tative sample. For a badly chosen initial proposal, this usually happens already at the first few 
iterations. By monitoring the file enc which is updated each iteration, an unsuccessful PMC run 
can be aborted. 
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The effective number of components can also be determined from any proposal file (mixjnvderis 



format) with the script neff .proposal . pi 



An additional diagnostic is the evolution of the proposal components with iteration. This il- 
lustrates whether the components spread out nicely across the high-posterior region and reach 
a more or less stationary behaviour, or whether they stay too concentrated at one point. The 
scripts proposaljTiean.pl (proposal_var.pl i read in the proposal information q' and plot 
the means (variances) as function of iteration t. 



3.3.2. Results 

PMC samples iter i/pmcsim 

This file contains the sample points. The first column is the (unnormalised) importance weight 
(log), the second column denotes the component number from which the corresponding point 
was sampled. Note that the Wciip points with highest weights are not considered in subsequent 
calculations (of moments, perplexity, evidence etc.). The next p columns are the p-dimensional 
parameter vector. Optionally, njed numbers of deduced parameters follow. 



Proposals iter_i/proposal 



The proposal used for the importance sampling in iteration / is in mixjtivdens format (Sect. A.3 1. 
The final proposal, updated from the sample of the last iteration, is proposal_f in. 



Mean and confidence intervals iter_i/mean 

This file contains mean and one-dimensional, left- and right-sided confidence levels (c.L). A 
c.l. of p% is calculated by integrating the one dimensional normalised marginal posterior starting 
from the mean in positive or negative direction, until a density of p%/2 is reached. PMC outputs 



c.l.'s for p = 63.27%, 95.45% and 99.73%. With the program cl_one_sided one-sided c.l.'s 
can be obtained. 



For post-processing, the program meanvar .sample outputs the same information (mean and 



C.l.) from an existing PMC sample, including possible deduced parameters 



Resampled PMC simulations iter_{niter-l}/sample 

If cosmo_pmc . pi has been run with the option -p, the directory of the final iteration contains the 
file of parameter vectors sample, which is resampled from the PMC simulation pmcsim, taking 
into account the importance weights. The resampled points all have unit weight. Resampling is 
a post-processing steop, it is performed by calling the R script sample_from_pmcsimu.R from 
cosmo_pmc . pi; this can also be done manually with any pmscim simulation. 
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Histograms iter_i/chi_j , iter_i/chi_j_k 

One- and two-dimensional histograms are written at each iteration to the text files chi_j and 
chi_j_k, respectively, where j and k, j < k, are parameter indices. Those histograms can be 
used to create Id- and 2d-marginals, using the script plot_contour2d.pl The bin number is 
set by the config entry nbinhist. 



In post-processing, use histograms_sample to produce histograms from a PMC sample. This 
can be useful if deduced parameters have been added to the sample. 



Covariance iter i/covar* . fin 

The parameter covariance and inverse covariance are printed to the files covar . fin and, respec- 
tively, covarinv . fin. The addition "+ded" in the file name indicates the inclusion of deduced 
parameters. The covariance matrices are in "mvdens" -format (see Sect.|A.3[). 



Evidence evidence 

This file contains the Bayesian evidence as a function of iteration. Before the first iteration, 
the Laplace approximation using the Fisher matrix is printed to evidence_fisher if the file 
fisher exists. At each iteration /, iter_i/evidence_covarinv contains the Laplace approx- 
imation of the evidence from the inverse covariance matrix of the sample iter_i/pmcsim. 



3.3.3. Deduced parameters 



Deduced parameters can be part of a PMC simulation. These parameters are not sampling pa- 
rameters, but they are deduced from the main parameters. For example, if Qm and Qa are 
sampling parameters of a non-flat model, the curvature Q.k - + can be a deduced param- 
eter. 

In most cases, deduced parameters are ignored while running CosmoPMC. They are usually 
added to the PMC simulation after the sampling, for example using a script. In the case of 



galaxy clustering, add_deduced_halomodel adds deduced parameters which depend on the 



sampling parameters but also on the underlying cosmology and halo model. 



A PMC simulation with deduced parameters added can be used as input to histograms_sample 



to create the histogram files, now including the deduced parameters. These can then in turn be 
read by and |plot_contour2d . pi to produce Id- and 2d-marginals, including the deduced pa- 
rameters. Alternatively, the PMC simulation with added parameters can be resampled using 



sample_f rom_pmcsimu . R from which plots can be created by plot_conf idence . R 
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3.3.4. Other files 

Maximum-posterior parameter max.logP 



max_post stores its estimate of the maximum posterior in this file. 



Fisher matrix fisher 



The final result of go_fishing the Fisher matrix in mvdens (Sect. A.3 1 format 



Log files log jnax_post , log_fish, log_pmc 



max.post go_fishing and cosmo_pmc each produce their corresponding log file. 



4. Cosmology 

The cosmology part of CosmoPMC is essentially the same as the stand-alone package nicae^ 
This excludes the external program camb and the WMAP likelihood library, which are called by 
CosmoPMC for CMB anisotropics. Further, CosmoPMC contains a wrapper layer to communi- 
cate between the PMC sampling and the cosmology modules. 



4.1. Basic calculations 

A number of routines to calculate cosmological quantities are included in the code. These are 

• Background cosmology: Hubble parameter, distances, geometry 

• Linear perturbations: growth factor, transfer function, cluster mass function, linear 3D 
power spectra 



Non-Unear evolution: fitting formulae for non-Unear power spectra (Peacock & Dodds 



T996| [Smith etaL]|2003] ), emulators Peitmann et"all|2009l [MTOt |Lawrence"eraL||20T0l ), 



halo model 

Galaxy clustering: HOD model 

Cosmic shear: convergence power spectrum, second-order correlation functions and de- 
rived second-order quantities, third-order aperture mass skewness 

CMB anisotropics via camb. 



' http : //www2 . iap . f r/users/kilbinge/nicaea 
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Table 1 : Extrapolation of the power spectra 



snonlinear 

linear 
pd96 

srai th® 3 , smi th© 3 _de 
coyotelO 



333.6 /zMpc-i 
333.6/1 Mpc-i 
333.6/1 Mpc-i 
2.416 Mpc"' 



"ext 

«s - 4 
-2.5 

Eq. (61), 

no extrapolation 



Smith et al. 



(2003) 



4.1.1. Density parameters 

Both the density parameters (Qx = Px/Pc) and the physical density parameters (a>x - ^xh^) are 
valid input parameters for sampling with PMC. Internally, the code uses non-physical density 
parameters (Ox). AH following rules hold equivalently for both classes of parameters. Note that 
physical and non-physical density parameters can not be mixed, e.g. Q.^ and ojk on input causes 
the program to abort. 

The parameter for massive neutrinos, Ovniass. is not contained in the matter density = + 

A parameter which is missing from the input list is assigned the default value, found in the 
corresponding cosmology parameter file (cosmo.par), unless there is an inconsistency with 
other input parameters. E.g., if Ode and Q.k are input parameters, Q^n is assigned the value 

OjYi — 1 — fide ~ ~ ^y,mass> 

to keep the curvature consistent with D.^. 
A flat Universe is assumed, unless (a) both and Qde> or (b) Q.^ are given as input parameter. 



4.1 .2. Matter power spectrum 



Usually, models of the non-linear power spectrum have a limited validity range in k and/or 
redshift. For small k, each model falls back to the linear power spectrum, which goes as Psik) oc 



k"\ For large k, the extrapolation as a power law Psik) oc ?iext is indicated in Table 4.1.2 
See for more details on the models. 



The Coyote emulator In the coyote 10 case, the power spectrum is zero for k > /cmax- The 
same is true for redshifts larger than the maximum of Zmax = 1- See Eifler (2011 1 for an alterna- 
tive approach. 

The Hubble constant h can not be treated as a free parameter. For a given cosmology, it has to be 
fixed to match the CMB first-peak constraint £a - ^dx^lr^ = 302.4, where t/is is the distance to 
last scattering, and is the sound horizon. This can be done with the function set Jl®_Coyote, 
see Demo/lensingdemo . c for an example. When doing sampling with non-physical density 
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parameters, h has to be set at each sample point. Alternatively, the physical density parameters 
can be sampled, where h is set internally to match the CMB peak. 

4.1.3. Likelihood 

Each cosmological probe has its own log-likelihood function. The log-likelihood function is 
called from a wrapping routine, which is the interface to the PMC sampler. In general, within 
this function the model vector is computed using the corresponding cosmology routine. The 
exception are the WMAP-modules where the Cf 's are calculated using camb and handed over to 
the log-likelihood function as input. 

4.2. Cosmic shear 

CosMoPMC implements second- and third-order weak lensing observables. 



4.2.1. Second-order 



The basic second-order quantities in real space for weak gravitational lensing are the two-point 
correlation functions ^± (2PCF) (e.g |Kaiser|1992[ ), 



1 r 



(6) 



Data corresponding to both functions (slensdata=xipin) as well as only one of them (xip, 
xim) can be used. The aperture-mass dispersion ( [Schneider et al.|1998| ) 



1 



(V) 



is supported for two filter functions Ugid-) = u('&/6)/6^ ( Schneider et al.|l998 Crittenden et al. 



polynomial (raap2poly): u{x) - -(1 - x ) 

n 



x^\ii{i-xy, 



Gaussian (map2gauss): u{x) - 
The top-hat shear dispersion ( |Kaiser|1992| ) 

1 

<lrl'>E,B(^) = ^ 

is used with slensdata = gsqr. 



In 



1 - 



e 2 . 



4Ji(^ 



(8) 
(9) 

(10) 
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Pure E-/B-mode separating functions ( [Schneider & Kilbinger 2007 1 are chosen with slensdata 
= decomp_eb. For the lack of analytical expressions for filter functions to obtain these real- 
space statistics from the convergence power spectrum, they are calculated by integrating over 
the 2PCF. The integral is performed over the finite angular interval [i^min; '^maxl- The prediction 
for the E-mode is 



E=- \ d&&[TA&)^+(&)±T-(&)^-m- 



(11) 



Two variants of filter functions are implemented: The 'optimized' E-/B-mode function Fu & 



Kilbinger (2010 1 for which the real-space filter functions are Chebyshev polynomials of the 
second kind, 



Y^a„Un{x); Unix) 



n=Q 



sin[{n + l)arccosx] 
sin(arccos x) 



(12) 



The coefficients a„ have been optimized with respect to signal-to-noise and the Qm-o"8 Fisher 
matrix. The function E is defined as a function of the lower angular limit i9'niin- The ratio rj of 
lower to upper limit, rj - j9'min/'9^max is fixed. 

The second variant are the so-called COSEBIs (Complete Orthogonal Sets of E-/B-mode Inte- 
grals; Schneider et al.|20T0 ). We implement their 'logarithmic' filter functions. 



rj^(^) 



In 



_d_ 

^mii 



/7+1 



Nn - rnj). 



(13) 



The coefficients Cnj are fixed by integral conditions that assure the E-/B-mode decomposition of 
the 2PCF on a finite angular integral. They are given by a linear system of equations, which 
is given in Schneider et al. ( [2010, ). To solve this system, a very high numerical accuracy 
is needed. The Mathematica notebook file $COSMOPMC/par_files/COSEBIs/cosebi .nb, 
adapted from Schneider et al. (20101, can be run to obtain the coefficients for a given j^min 
and 7?max- An output text file is created with the zeros r,„- and amplitudes N^- The file name is 
cosebi_tplog_rN_[Nmax] .[thmin] _[thmax] , where Nmax is the number of COSEBI modes, 
thmin and thraax are the minimum and maximum angular scale ^^mn and i^max. respectively. 
For a given ^y^m and i?max> specified with the config entries thjnin and thjnax, CosmoPMC 
reads the corresponding text file from a directory that is specified by path. A sample of files 
with various scales are provided in $COSMOPMC/par_files/COSEBIs. 

The COSEBIs are discrete numbers, they are specified by an integer mode number n. 

In both cases of pure E-/B-mode separating statistics, the function r_ is calculated from T+ 
according to |Schneider et al.| ( |2002| ). 

The additional flag decomp_eb_filter decides between different filter functions: 
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decomp_eb_filter Reference 



Filter function typ 



FKIQ.SN 

FK10_FoM_etal® 
FK10_FoM_eta5® 
COSEBIs_log 



Fu & Kilbinger 


(2010 


Fu & Kilbinger 


(2010 


Fu & Kilbinger 


(2010 


Schneider et al. 


(2010 



optimized Signal-to-noise 1 /50 

optimized Fisher matrix 1/10 

optimized Fisher matrix 1 /50 
logarithmic 



The convergence power spectrum with covariance matrix can be used with the flag slensdata 
= pkappa. 



4.2.2. Third-order 

We implement the aperture-mass skewness ( |Pen et al.|2003t [Jarvis et al.|2004| [Schneider et al. 



2005| l with the Gaussian filter (eq. ??). There are two cases: 
slensdata = mapSgauss 



The 'generalised' skewness (^ap)( > ^3) ( Schneider et al. 
slensdata 



2005 1 with three filter scales. 



map3gauss_diag 
The 'diagonal' skewness ^M^p)^ 6) using a single aperture filter scale. 

TODO: equations 



4.2.3. Second- plus third-order 

A joint data vector of second- and third-order observables can be used in CosmoPMC. The 
covariance is interpreted as a joint block matrix, with the second-order and third-order auto- 
covariances on the diagonal, and the cross-correlation on the off"-diagonal blocks. The possible 
scenarios are: 

• slensdata = map2gauss_map3gauss 

Gaussian aperture-mass dispersion and generalised skewness. 

• slensdata = map2gauss_map3gauss_diag 
Gaussian aperture-mass dispersion and diagonal skewness. 

• slensdata = decomp ebjnap3gauss 

Log-COSEBls and generalised aperture-mass skewness. The flag decomp_eb_f ilter has 
to be set to COSEBIs_log. 

• slensdata = decomp_eb_map3gauss_diag 

Log-COSEBls and diagonal aperture-mass skewness. The flag decomp _eb_f ilter has to 
be set to COSEBIs_log. 
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The first two cases use the same fiher for second- and third-order, and provide therefore a con- 
sistent measure for both orders. The last two cases use the optimal E-/B-mode function known 
for second order. 



4.2.4. Covariance 



The covariance matrix is read from a file, and the inverse is calculated in CosmoPMC. The matrix 
has to be positive definite. An Anderson-Hartlap debiasing factor is multiplied to the inverse 



(Anderson 2003} Hartlap et al. 2007 1, which is specified with the config entry corr_invcov.This 
can also be used to rescale the covariance, e.g. to take into account a different survey area. Set 
this value to unity if no correction is desired. 

The covariance is either taken to be constant and not dependent on cosmology. In that case. 



set scov_scaling to cov_const. Or the approximated schemes from Eifler et al. (2009 1 are 
adopted, see ? for the implementation. In that scheme, the shot-noise term D is constant, the 
mixed term M is modulated with Om and crs using fitting formluae, and the cosmic-variance 
term V is proportional to the square of the shear correlation function. This scheme is available 
for slensdata = xipm. The three covariance terms have to be read individually. The entry 
covname, which for scov_scaling = cov_const corresponds to the total covariance matrix, 
now specified the file name of cosmic-variance term, covnameJl the name of the mixed term, 
and covname_D the name of the shot-noise term. 



4.2.5. Reduced shear 

The fact that not the shear y but the reduced shear g = y/{l-K) is observable leads to corrections 
to the shear power spectrum of a few percent, mainly on small scales. These corrections are 



either ignored, or modelled to first order according to Kilbinger (2010 1. This is controlled in 
the lensing parameter file (cosmo.lens.par). The parameter range where the reduced-shear 
corrections are valid are indicated in Tabled 

4.2.6. Angular scales 

The flag s format describes the mapping of angular scales (given in the data file) and 'effective' 
scales, where the model predictions of the shear functions are evaluated: 

1. sformat = angle_center: The effective scale is the same as given in the data file, 

2. sformat = angle jnean: The model is averaged over a range of scales [%, 9i] given in 
the data file. 
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Table 2: Parameter limits where the reduced- shear corrections are valid (from Kilbinger 20101 



a 


Parameter 


lower 


upper 


1 




0.22 


0.35 


2 




0.33 


1.03 


3 


W 


-1.6 


-0.6 


4 




0.005 


0.085 


5 




0.61 


1.11 


6 


0"8 


0.65 


0.93 


7 


n,s 


0.86 


1.16 



3. sformat = angle_wlinear: The model is the weighted average over a range of scales 
[00, 9i], where the weight is w = 0/arcmin. 

4. sformat = angle.wquadr: The model is the weighted average over a range of scales 
[6o, 9i], where the weight is w = fl'i(0/arcmin) + a2i9 / arcmin)^ . 

The first mode (angle_center) should be used for aperture-mass, shear rms and 'ring' statistics, 
since those quantities are not binned, but instead are integrals up to some angular scale 6. For the 
correlation functions, in particular for wide angular bins, one of the last three modes is preferred. 
The quadratic weighting (angle.wquadr) corresponds to a weighting of the correlation function 
by the number of pair^ This mode was used in the COSMOS analysis ( Schrabback et al. 2010 1. 



4.3. SNIa 

The standard distance modulus (schi2mode = chi2_simple) for a supernova with index / is 

f^B,i = mi l - ^ + «(^; - 1) -^Ci. (14) 

where the quantities measured from the light-curve fit are the rest-frame B-band magnitude ., 
the shape or stretch parameter si, and the color c,-. The universal absolute SNIa magnitude is M, 
the linear response parameters to stretch and color are a and /?, respectively. The -function is 

where di^ is the luminosity distance and Zi the redshift of object /. The contributions to the 
total error for object / are: (1) The light-curve parameter variance (T^(jubj) = Q^^^Q^ the 



'p. Simon, private communication 
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parameter vector 02 = {l,a,/3) and the covariance W2 of the data vector {m*g., Si,Ci). (2) The 
peculiar velocity uncertainty o"pv,, = 5/ In 10 • Vp/(cz,)- (3) The intrinsic absolute magnitude 
scatter o-{n\.- 

The Hubble parameter is absorbed into the absolute magnitude which we define as M = M - 

51ogio/l70- 



The form of this log-likelihood function has been used in Astier et al. (2006 1. 

The following variations of the distance modulus and log-likelihood are implemented: 

• schi2mode = chi 2 _Thetal: The is extended to include photometric zero-point un- 



certainties, see Kilbinger et al. (20091. 

schiZmode = chi2_ThGta2_denora_fixGd: The parameters a and yS in the denominator 



of ( 15 1 are fixed and kept constant during the Monte-Carlo sampling. 

schi2modG = chi2 jio.sc: The stretch and color parameters are ignored, the distance 
modulus is hbj = m*^ - M. 

schi2mode = chi2_betaz: Instead of a single parameter, the color response is redshift- 
dependent, /3 ^ j3 + p^zi- 

chi2 dust: Intergalactic dust absorption is taken into account in the distance modulus, 
see 



Menard etal. (2010 1. 



The covariance matrix W2 of the data vector (m*^ ., Si, c,) depends on the parameters a and p. In 
a Bayesian framework, this leads to an additional term j log det W2 in the log-likelihood func- 
tion. Taking into account this parameter-dependent term leads however to a biased maximum- 
likelihood estimator, in particular for a and Therefore, it is recommended to not include this 
term. Use the flag add_logdetCov = 0/1 in the configuration file to disable/enable this term. 



4.4. CMB anisotropies 

The full CMB ani sotropics are handled externally: The Q's are calculated by calling camtj^ 
( [Lewis et aLl|2000| , the W MAP likelihood functi on (3'^''-, 5*- and 7*-year) is computed using 
the WMAP public code|^ ( Dunkley et al. 2009 1. The maximum { up to which the C('s are 
calculated and used in the likelihood can be determined in the configuration file. An ^^ax - 2000 
is recommended for high precision calculations. 

The power spectrum from the Sunyaev-Zel'dovich (SZ) effect can be added to the Ce's, multi- 
plied with an amplitude A as free parameter. The predicted SZ power spectrum is taken from 
Komatsu & Seljak (20()2|l. This model has been used in the 3-, 5- and 7-year analyses of the 
WMAP data ( |Komatsu et al.|2011| l. 



''J. Guy, private communication 
' http : //camb . info 



^http : //lambda ■ gsf c . nasa . gov 
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Alternatively, the WMAP distance priors ( Komatsu et al.|2009 1 can be employed. 



4.5. Galaxy clustering 
4.5.1. Halomodel and HOD 



The theoretical model of galaxy clustering is the one used in |Coupon et al. (2012 ); see this paper 
for details of the model and further references. 



As the basis to describe galaxy clustering, we implement the halo-model as reviewed in ( |Cooray| 
& Sheth||2002 1, which accounts for the clustering of dark-matter halos. On top of that, a halo 
occupation distribution (HOD) function ( jBerlind & Weinberg|2002 Kravtsov et al. 2004[ Zheng 



et al. 2005;) is the prescription of how galaxies populate those halos. This function is the number 
of galaxies N in a. halo of mass M. With the flag hod = berwein02_excl, this number is 
expressed as the sum of central (A^c) plus satellite (A^s) galaxies. 



N{M) - A^c(M) x[l+ N,{M)] 



with 



ndM) 
n,{M) 



1 -I- erf 



logio M - logio M„ 



O'logM 

(^f; ifM>Mo 







else 



(16) 

(17) 
(18) 



We further compute the galaxy two-point correlation function ^(r) and its angular projection 
w{9) using the redshift distribution provided by the user, as well as the galaxy number density 
(for a full description of the model see Coupon et al.|2012 1. To prevent haloes from overlapping, 
we implement the halo exclusion formalism as described in Tinker et al.| ( |2005| ). 

For the halo bias, three options are available: 
• shale bias = bias sc 



Bias expansion from the spherical collapse model, see e.g. eq. (68) from Cooray & Sheth 



(20021. 



shalo_bias = bias.tinkerOS 



Bias calibrated with numerical simulations. Tinker et al. (2005 1 eq. (Al). 
shalo bias = bias tinkerlO 



Updated bias fitting formua from Tinker et al. (2010 1, eq. (6) and Table 2. 
The mass function describes the number of halos for a given mass and redshift. It is defined as 

d« Po v/(v) dy 



dlnM M 



dlnM' 



(19) 
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where v(M, z) - dc{z)/[D+{z)o-{M)] is a measure of the overdensity with (t(M) being the rms 
matter fluctuation in a top-hat window containing the mass M. Pq = QmPcfl is the mean density 
of matter at the present day. 

The following mass functions are implemented, via the flag smassfct: 
• From the spherical/eliptical collapse model: 



- ps: p = 0,q^l ( [Press & Schechter|1974l ) 

- St: p = 0.3, <7 = 0.75 ( |Sheth & Tormen|1999 1 

- st2: p = 0.3, <7 - 0.707 ( |Sheth & Tormen|1999| l 
From numerical simulations: 



y/(v) -/(0-) = 0.315 exp[-|ln(o--i +0.61p-^] 
jSl: ( [Jenkins et al.|2001| ) 



(20) 



(21) 



The dark-matter halos have the density profile 

p(r)=p,[{r/r,m+r/r,)'-"l\ (22) 

For slopes unequal to the [Navarro et al.[ ( [l997[ ) value of a = 1, closed expressions for the Fourier 
transform of p do not exist, and the code will be slower. 

The concentration parameter is given by 

" M 



c{M,z) = 



l+z 



(23) 



following Takada & Jain (2003 1. The parameters cq and yS can be chosen freely in the halomodel 
parameter file halomodel . par. 

The log-likelihood function is the sum of the contribution from the angular correlation function 
and the galaxy number density Wgai: 

Ij^obs _ ^modell'^ 



ngal 



where n^'j"^^' is estimated at the mean redshift of the sample. 



The number of galaxies (second term in eq. 24 1 can be included in the following way, with the 
config flag sngal_fit_type: 
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• ngal_lin_fit: linear (standard; according to the above equation) 

• ngal_log_fit: logarithmical 

• ngal Jio_fit: no inclusion, second term is omitted 

• ngal_lin_fit_only: exclusive, first term is omitted 

4.5.2. Deduced parameters 

The following deduced parameters can be computed: 

• Mean galaxy bias 



(z)^ JdMbi,{M,z)n{M,zy-^, (25) 



«sai(z) ' 

where is the halo bias, and 



-I 



«gai(z)= N(M)n(M,z)dM (26) 



is the total number of galaxies. 
Mean halo mass 



N(M) 

{M^^o){z)= \ dMMn{M,z)^-{. (27) 

«gal(z) 



Fraction of satellite galaxies 



/s(z)- l-/c(z); /c(z)- JdMn(M,z)^^^. (28) 



N,{M) 

«gal(z) ' 



Use the program add_deduced Jialomodel to add those deduced parameters to a PMC sample. 



See the example config file conf ig_pmc_ded in Demo/MCJ)emo/H0D/CFHTLS-T®6. 
4.5.3. Clustering data 

The angular two-point correlation function w{0) is implemented, with the flag shalodata = 
wof theta. The measured (input) data Wmes is corrected for the integral constraint, via 

w{9) - wra^m + wc, (29) 

assuming that the measured correlation function can be fit by a power law 

Wmes(e)«A^(r'5-C). (30) 



The program haloplot outputs the correlation functions w{6) and ^(r), the HOD function N{M), 



and deduced parameters for given HOD input parameters. 
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4.5.4. Comoving volume 



The comoving volume is needed to calculate the comoving number density of galaxies, follow- 
ing from the halomodel and the HOD parameters. There are two possibilities to calculate the 
comoving volume Vq- First, if Zmin and z^ax are larger than zero in the HOD parameter file 
halomodel .par (see Table |7]l, Vq is computed between those two redshifts. Second, if both 
numbers are < 0, Vq is weighted by the redshift distribution n(z), see e.g. eq. (28) in Ross & 



Brunner (2009 1. In this weighting, the maximum value of n{z) is set to unity. 



4.6. BAO 



BAO constraints are implemented with two distance measures: 
• smethod = di stance 



The distance parameter A is defined in Eisenstein et al. (2005 1 as 



A(z) 



c/Ho z 



where 



Dwiz) 



flMz)] 



cz 



1/3 



(31) 



(32) 



is the spherically averaged distance to redshift z- 
smethod = distance_d_z 

The distance parameter d is the ratio of sound horizon at drag epoch Zd to spherically 
averaged distance (e.g. |Percival et al.|2007[ ). 



d{z) 



rsiZd) 
DvizY 



(33) 



We use the fitting formula for the drag redshift Zd from Eisenstein & Hu ( |1998 1 and cal- 
culate the sound horizon as the distance a sound wave can travel prior to Zd by numerical 
integration. 



4.7. Redshift distribution 

Some of the cosmology modules require a redshift distribution, for example lensing and HOD. 
Table[3]lists the implemented redshift distributions n{z), via the flag nofz. 

Each redshift bin can have a different type. The syntax for a redshift bin file is described in 
Appendix |A.1.5 
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Table 3: Redshift distribution types 



nofz Description 



parameter list 



hist 

single 

ludo 

jonben 

yinink 



Histogram I^ta ' Mz,;z.,^ 

Single redshift 6b{z - Zo) 
Fitting function {z/zoT exp [- (z/zo)^] 
zV{z 



(z" + z"^)/(. 



z' + c) 



(see text) 

Z0,Z() 

Zmirn Zmax> ^^tf^t Zo 
Zmiri) Zmax> b, C 
Zmin> Zmax> ^J, b, C 



All redshift distributions are internally normalised as 



/^^max 

Zmin 



dzn{z) - 1. 



(34) 



4.8. CMB and the power spectrum normalisation parameter 

The power spectrum normalisation parameter taken as input for camb is A^, which is the ampli- 
tude of curvature perturbations at the pivot scale ko = 0.002 Mpc~^ For lower-redshift probes 
such as lensing or HOD, the normalisation is described by erg, the rms fluctuation of matter in 
spheres of 8 Mpc/h. To combine those probes in a PMC run, has to be an input parameter, 
and (T8 a deduced parameter. CMB has to come first in the list of data sets so that camb can 
calculate erg, which in turn is handed over to the lensing likelihood. 



4.9. Parameter files 

Tables |4]-|6]list the contents of the parameter files for basic cosmology, lensing, SNIa and HOD. 
Proto-types can be found in $COSMOPMC/par_files. These files specify the default values of 
parameters and flags. These default values are over-written if any of those parameter is used for 
Monte-Carlo sampling. 



5. The configuration file 

The programs max_post, go_fishing, cosmo_pmc, and cosmo_mcmc read a configuration file 
on startup. Each configuration file consist of two parts: 

The first, basic part is common to all four config file types (Table [9]l. It consists of (1) the 
parameter section, (2) the data section and (3) the prior section. The data-specific entries in the 
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Table 4: Basic cosmology parameter file (cosmo . par) 



Omegajn Q.^ 

Omega.de Qde 

w0_de Wo 

wl_de wi 

h_100 h 

Oraega_b Qb 

Oraegajiujnass Q. 
N_eff Jiujnass 



n_spec 
snonlinear 



stransfer 



sgrowth 



sde param 



norinmode 
ajTiin 



v,mass 
v.mass 



normalization erg 



linear 
pd96 
smith® 3 
smith®3_de 



Matter density, cold dark matter + baryons 
Dark-energy density (if w = - 1 , corresponds to Qa 
Dark-energy equation-of-state parameter (constant term) 
Dark-energy equation-of-state parameter (linear term, 
see sde_param) 

Dimensionless Hubble parameter 
Baryon density 

Massive-neutrino density (so far only for CMB) 
Effective number of massive neutrinos (so far only for 
CMB) 

Power- spectrum normalisation at small scales (for 
normmode==0, see below) 
Scalar power-spectrum index 
Power spectrum prescription 
Linear power spectrum 



Peacock & Dodds 



Smith et al. 



(20031 



(19961 



Smith et al. (2003 1 + dark-energy correction from 



1 cosmo . org 



coyotel® 'Coyote Universe', Heitmann et al. (20091, Heitmann 



et al. (20101, Lawrence et al. (20101 



Transfer function 



bbks 

eisenhu 

eisenhu_osc 

heath 
growth.de 



jassal 
linder 



Bardeen et al. 



(19861 



Eisenstein & Hu 



Eisenstein & Hu 



(1998 



(1998 



'shape fit' 

with BAO wiggles 



Linear growth factor 

(11977') fitting formula 



Heath 



Numerical integration of differential equation for 6 (rec- 
ommended) 

Dark-energy parameterisation 

w(a) = Wo + wia{\ - a) 

w(a) = Wo -I- wi(l - a) 

Normalization mode. 0: normalization^crg 

Minimum scale factor 



data section are listed in Table [TT] 



The second part is type-specific. See Table 10 for the PMC part, and Table 13 for the MCMC 
part. Example files can be found in subdirectories of $COSMOPMC/Demo/MC_DEMO. 
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Table 5: Weak lensing parameter file (cosmo_lens . par) 



cosmo_file 
nofz.file 
redshift module" 
stomo 



sreduced 



q_mag_size 



tomo_all 

tomo_auto_only 

tomo_cross_only 

none 

K10 

1 



Basic cosmology file name (cosmo . par) 

Redshift distribution master file 

(see Table [sjl 

Tomography correlations 

All correlations 

Only auto-correlations (//) 

Only cross-correlations (/ j) 

Reduced-shear treatment 

No correction 



Fitting-formulae from Kilbinger (20101 



Magnification-bias coefficient, q 
(see,Kilbinger|2010[ eg. 16) 



Table 6: SNIa parameter file (cosmo.SN . par) 

cosmo.f ile Basic cosmology file name (cosmo . par) 

Theta2 -Ma - P Pz Distance modulus parameters 



To create a config file of type max.post or go_fishing from a PMC config file, the script 
|conf ig_pmc_to_max_and_f ish . pl| can be used. 

Some flags are handled internally as integers (enumerations), but identified and set in the config 
file with strings. The corresponding key word carries the same name as the internal variable, 
preceded with an 's', e.g. the integer/string pair lensdata/slensdata. 

The prior file, indicated if desired with the flag sprior, is a file in mvdens format. It specifies 
a Gaussian prior with mean and covariance as given in the file. Note that the covariance and not 
the inverse covariance is expected in the file. 
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Table 7: HOD parameter file (halomodel . par) 



cosmo.file 

nofz_file 

redshift module" 

zrain 

zmax 

alpha_NFW 

C0 

beta_NFW 
smassfct 



Mjtiin 
Ml 

m 

sigma_logjyi 

alpha 

shod 



Zmin 
Zmax 

a 
CO 



ps 

St 

st2 
j01 

^min 

Ml 
Mo 

O'logM 

a 

berwein®2 Jiexcl 



Basic cosmology file name (cosmo . par) 
Redshift distribution master file 
(see Table [8| 

Minimum redshift (-1 if read from nzf ile) 
Maximum redshift (-1 if read from nzf ile) 
Halo density profile slope (a = 1 for NFW) 
Concentration parameter at z = 
Concentration parameter slope of mass depen- 
dence 

Halo mass function type 

(|Press & Schechter|l974|), p = 0,q=\ 



( Sheth & Tormen 



( Sheth & Tormen 



1999 



1999 



,p = 03, q = 0.75 

, p = 0.3, q - 0.707 

( [Jenkins et al.|2001[ ) 
Minimal mass for central galaxies [/j~'Mq] 
Scale mass for satellites [H^^Mq] 
Minimum mass for satellites [H'^Mq] 
Logarithmic dispersion for central galaxies 
Slope for satellite mass dependence 
HOD type 



Berhnd & Weinberg (2002 1 with halo exclusion 



Table 8: Redshift module file (nof z . par) 



Nzbin 

snzmode nz_read_from_files 
nzf ile fi, f2, fNzbin 



Number of redshift bins 
File mode 

File names. See Appendix 



A.1.5 



for the file syntax. 
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5. The configuration file 



Table 9: Basic, common part of the configuration file 

version double Config file version. Upwards compatibility (config file 

version > CosmoPMC version) cannot be guaranteed. 
Downwards compatibility (config file version < Cos- 
moPMC version) is most likely ensured. 



Parameter section 



npar 
n ded 



spar 



min 
max 



mteger 
integer 



string 



npar+n_ded doubles 
npar+n_ded doubles 



Number of parameters 

Number of deduced parameters. The deduced parame- 
ters are not sampled but deduced from the other param- 
eters and written to the output files as well 
Parameterisation type, necessary for the wrapping into 
the individual posterior parameters and for plotting, see 
Table 



12 for possible parameters 



Parameter minima 
Parameter maxima 



ndata integer 
sdata string 

sdata string 



Data section 

Number of data sets 
Data set 1 

Data set ndata 



Prior section 



sprior 


string 


Prior file name (' 


'-" for no prior) 


[nprior 


integer 


If sprior t 


Number of parameters to which prior 






applies] 




[indprior 


npar x {0, 1) 


If sprior i= 


Indicator flags for prior parameters] 
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5. The conGguration Hie 



Table 10: PMC part of the configuration file 



nsample 


integer 


Sample size per iteration 


niter 


integer 


Number of iterations 


fsfinal 


integer 


Sample size of final iteration is fsfinal X nsample 


niter 


integer 


Number of iterations (importance runs) 


nclipw 


integer 


The nclipw points with the largest weights are discarded 


Proposal section 


d£ 


double 


Degrees of freedom (df=-l is Gaussian, df=3 is 'typical' 






Student-t) 


ncomp 


integer 


Number of components 


sdead corap 


string 


One of 'bury' , 'revive' 


sinitial 


string 


Proposal type (one of fisher_rshift, fisher_eigen, 






file, random_position) 


f Shi ft'' 


double 


Random shift from ML point ~ U{-r, r); 






r - fshift/Cmax-min) 


£var° 


double 


Random multiplier of Fisher matrix 


prop_ini_name* 


string 


File name of initial proposal 


fmin*^ 


double 


Components have variance ~ f7(fitiin, (max - min)/2) 


Histogram section 


nbinhist 


integer 


Number of density histogram bins 
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5. The configuration file 



Table 11: Data-specific entries in the configuration file's data section 



Weak gravitational lensing Lensing 



slensdata 



th_min* 
thjnax'' 
path^ 
s format 



ar 

datname 
scov.scaling 
covname 
covname.M'' 
covname _D'' 
corr invcov 

Nexclude 

exclude^ 



string 



sdGcomp_eb_filter'' string 



double 
double 
double 
string 



double 
double 

string 
string 
string 
string 
string 
double 

integer 



xip, 



Nexclude integers 



Data type, one of xipm, 
xim, map2poly, map2gauss, 
gsqr, decomp eb, pkappa, 
mapSgauss, map3gauss_diag , 
map2gaussjnap3gauss , 
map2gauss_map3gauss_diag , 
decomp_ebjnap3gauss , 
decomp_eb_map3gauss_diag 
One of FK10_SN, FK10_FoM_etal0 , 
FK10_FoM_eta5Q, COSEBIs_log 
Minimum angular scale 
Maximum angular scale 
Path to COSEBIs files 

Data format of angular scales, one 
of angle_center , anglejnean, 

angle_wliriear , angle_wquadr 
Linear weight 

Quadratic weight, w = al • 0/arcmin + a2 • 

(^/arcmin)^ 

Data file name 

One of cov.const , cov_ESHQ9 
Covariance file name 
Covariance mixed term file name 
Covariance shot-noise term file name 
Correction factor for inverse covariance ML 



estimate, see jHartlap et al. (2007 1 

Number of redshift bin pairs to be excluded 

from analysis 

Indices of redshift pairs to be excluded 



Parameter file name, e.g. cosmo_lens 
Additional prior, one of none (recom- 
mended), unity , de_conservative 



model file 
sspecial 



strmg 
string 
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5. The configuration file 



Supernovae type la SNIa 



datname 

datformat 

schi2mode 



ThetaZ.denom'' 
zAVjiame^ 
datname_beta_d'' 
add logdetCov 



model-file 
sspecial 



string 
string 
string 



2 doubles 
string 
string 
integer 

string 
string 



Data file name 

Data format, SNLS_f irstyear 

X'^ and distance modulus estimator type (one 
of chi2_simple , chi2_Theta2_denora_fixed, 

chi2 betaz, chi2 dust, chi2 residual) 

Fixed a,/3 in ;^f^-denominator 
File with Ayiz) table 

Prior file (mvdens format) on ySd ("-" if none) 

1 if 0.5 log det Cov is to be added to log-likelihood, if not 



(recommended; see Sect. 4.3 1 



Parameter file name, e.g. cosmo_SN 

Additional prior, one of none (recommended), unity, 
de_conservative 



Table 11: Data-specific entries in the configuration file's data section (continued). 



CMB anisotropics CMB 



scamb_path 
data.path 



Cl_SZ_file 
Imax 

accurate 
model file 
sspecial 



strmg 
string 



string 

integer 

on 

string 

string 



/path/to/ scamb 

/path/to/wmap-data. This path should contain the direc- 
tory data with subdirectories healpix_data, highl, lowlP, 
lowlP 

File with SZ correction angular power spectrum ("-" if none) 
Maximum £ for angular power spectrum 
Accurate reionisation and polarisation calculations in camb 
Parameter file name, e.g. cosmoDP .par 

Additional prior, one of none (recommended), unity, 
de_conservative 



WMAP distance priors CMBDistPrior 

datname string Data (ML point and inverse covariance) file 
model.f ile string Parameter file name, e.g. cosmo_lens . par 

sspecial string Additional prior, one of none (recommended), unity, 
de_conservative 
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5. The configuration file 



Galaxy clustering (HOD) GalCorr 



shalodata 
shalomode 
datname 
covname" 
corr invcov 

delta 

intconst 

area 

sngal_fit_type 



ngal* 
ngalerr^ 
model.file 
sspecial 



string 
string 
string 
string 
double 

double 
double 
double 
string 



double 
double 
string 
string 



galcorr_log 



Data type, wof theta 

type, one of galcorr_var , galcorr_cov , 
Data (+variance) file name 
Covariance file name 

Correction factor for inverse covariance ML estimate, see Hart 



lapetal. (20071 



Power-law slope 5, for integral constraint 
Integral constant C 
Area [deg^] 

Likelihood type, inclusion of galaxy number. One 

of ngal_lin_fit , ngal_log_fit , ngal_no_fit, 

ngal_lin_fit_only 

Number of observed galaxies 

Error on the number of observed galaxies 

Parameter file name 

Not used for HOD, set to none 



Table 11: Data-specific entries in the configuration file's data section (continued). 



Baryonic acoustic oscillations BAO 



smethod string BAO method, one of distance_A, distance_d_z 
datname string Data -i- covariance file name (mvdens format) 
model.file string Parameter file name, e.g. cosmoDP .par 
sspecial string Additional prior, one of none (recommended), 
de_conservative 



unity, 



Table 12 contains a list of input parameters, which can be given as strings to the spar key in the 
config file. 
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5. The configuration file 



Table 12: Input parameters 
Name Symbol Description 



Basic cosmology (some of them given in cosmo . par) 

Omegajn Q.^ Matter density, cold dark matter + baryons 

omegajti a»m 

OmGga_b Q.], Baryon density 

omega.b Wb 

lQO_omGga_b 100 x aj\, 

Omega.de Q^c Dark-energy density (if w = - 1 , corresponds to 

omega_de a»de 

Omegajiujnass f^y,mass Massive-neutrino density (so far only for CMB) 

omegajiujnass <i»v,mass 

Omega.c Q.c Cold dark matter 

omega_c a»c 

OmegaJC Q.k Curvature density parameter 

omegaJC ojk 

wO_de Wo Dark-energy equation-of-state parameter (constant term) 

wl_de wi Dark-energy equation-of-state parameter (linear term, see 

sde_param) 

h_100 h Dimensionless Hubble parameter 

N_ef f JiujTiass A'^eff.v.mass Effective number of massive neutrinos (so far only for CMB) 

sigma_8 erg Power-spectrum normalisation at small scales 

Delta_2_R Power-spectrum normalization at large scales (CMB) 

n_spec ?is Scalar power-spectrum index 

alpha_s as Running spectral index (so far only for CMB) 

n_t Hi Tensor power-spectrum index 

r r Tensor to scalar ratio 

ln_r In r 

tau T Optical depth for reionisation 

A_SZ Asz SZ-power spectrum amplitude 
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6. Post-processing and auxiliary programs 



Table 12: Input parameters (continued) 



SNIa-specific (some of them given in cosmo_SN . par) 

M M - log 10 hjQ Universal SNIa magnitude 

alpha a Linear response factor to stretch 

beta j3 Linear response factor to color 

beta_z ySj. Redshift-dependent linear response to color 

beta_d j3d Linear response to the color component due to intergalactic dust 



Galaxy-clustering- 


specific (some of them 


given in halomodel . par) 


MjTiin 




Minimum halo mass for central galaxies 








loglOJljnin 


logioMmi„/(Mo/i-i) 




M_l 


Ml 


Scale mass for satellite galaxies [Mq/j"^] 


logl®Jl_l 


logio[Mi/(Mor^)] 




M_® 


Mo 


Minimum halo mass for satellite galaxies 






[Mo/j-i] 


logl®Jl_0 


logioMo/(Mori) 




sigma_log_M 


ClogM 


Dispersion for central galaxies 


alphaJialo 




Slope of satellite occupation distribution 


MJialo.av* 


<Mh) 


Average halo mass [Mq/j"^] 


logl®JlJialo_av 


* logio{M^I{MQh-')) 




b_halo_av* 


(bh) 


Average halo bias 


N_gal_av* 




Average galaxy number per halo 


fr.sat* 


/s 


Fraction of satellite galaxies to total 


ngal_den* 


rig 


Comoving galaxy number density [Mpc'^h^] 


logl®ngal_den* 


logio?lg 





6. Post-processing and auxiliary programs 

All scripts described in this section are located in $COSMOPMC/bin. 
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6. Post-processing and auxiliary programs 



6.1. Plotting and nice printing 
6.1.1. Posterior marginal plots 

Marginals in Id and 2d can be plotted in two ways, using (1) plot_contour2d.pl or (2) 
plot_confidence .R. The first is a perl script calling yorick for plotting, the second is an 
R script. The second option produces nicer plots in general, in particular, smoothing workes 
better without producing over-smoothed contours. Further, filled contours with more than one 
data set aie only possible with the R option, yorick can only combine several plots with empty 
contours. The computation time of the R script is however much longer. 



1. plot_contour2d.pl creates Id and 2d marginals of the posterior, from the histogram 
files chi2_j and chi2_j_k. 

To smooth Id and 2d posteriors with a Gaussian, use plot_contour2d.pl -n -g 
FACTOR. The width of the Gaussian is equal to the box size divided by FACTOR. It is rec- 
ommended to test the smoothing width FACTOR by setting it to a negative number which 
causes both smoothed and unsmoothed curves being plotted. This can reveal cases of 
over-smoothing. If contours have very different width in different dimension, the addition 
option -C uses the PMC sample covariance (from the file covar+ded . fin) as the covari- 
ance for the Gaussian. For the final plot, replace -FACTOR with FACTOR to remove the 
unsmoothed curves. Remove the option -n to add color shades to the 2d contours. 

The file log_plot contains the last plot command with all options. This can be used to 
reproduce and modify a plot which has been generated automatically by other scripts, e.g. 
[cosmo pmc .pi] 



2. plot_conf idence . R creates Id and 2d marginals of the posterior, from the re-sample file 
sample. 

Smoothing is done with a kernel density estimation using the R function kde2d. The 
kernel width can be set with the option -g. The number of grid points, relevant both for 
smoothing and filled contours, is set with -N. Use both -i and - j options to only plot the 
2D marginals of parameters iand] to save computation time. 

6.2. Mean and confidence intervals 

From a "mean" output file, containing parameter means and confidence levels, one can create a 



ps/pdf file using the command mean2eps . pi 



This is equivalent to the following steps (see also essential_cosmo_pmc_run.pl l: 



|meanvar2tab . pi creates a table with parameter names and values formatted in TgX- 



format. 



tab2tex . pi wraps a ETgX table header around the table. 
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6. Post-processing and auxiliary programs 



• txt2tex.pl wraps a I^TgX header around the file. 
Example: 

> meanvar2tab.pl -s 1 -p 2 -e iter_9/mean > mean. tab 

> tab2tex.pl -s 1.25 mean. tab > mean in.tex 

> txt2tex.pl mean_in.tex > mean.tex 



6.2.1. PMC proposal 



proposaljnean . pi (proposal_var . pi I creates plots of the proposal component's means (vari- 



ances) as function of the iteration. 



6.3. Importance sampling 

A PMC simulation file (pmcsim) from an earlier PMC run, corresponding to a sample from 
posterior pi, can be used to do importance sampling with another posterior p2. For that, simply 
replace the data section of the earlier config file with the corresponding data section of posterior 



P2. The command importance_sample creates a new PMC simulation which corresponds to a 



sample under the posterior product p\ ■ P2- 
6.4. Bayesian evidence, Bayes' factor 



evidence . pi calculates and prints the evidence from a PMC simulation file. The same infor- 



mation is printed to the file evidence during a PMC run. 



bayes.f actor . pi prints Bayes' factor between two PMC runs together with the Jeffrey scale. 



evidence_list . pi prints a list of evidences for a number of PMC runs. 



6.5. Reparameterisation 



remap . sh swaps and removes parameters from a MCMC or PMC run. The histogram files, 
mean and covariances are remapped. This is useful if different runs are to be reduced to a 
common parameter set for comparison or joint plotting. The removal of parameters is equivalent 
to marginalisation over the corresponding parameter subspace. 

For example, suppose there is a SNIa run in directory Sn, and a lensing run in Lensing. SNIa 
has the following parameters: 

Omegam Omegade w0de M alpha beta 
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6. Post-processing and auxiliary programs 



Lensing has the parameters: 

Omegam sigmaS wOde Omegade hlO® 

In Sn, create the file remap . dat with the hne 

® 1 2 

In Lensing, create the file remap . dat with the line 
3 2 

In both directories run the command 

> remap. sh -i iter_<niter-l> 

which creates sub-directories remap containing symbolic links and/or copies of histogram files 
to/from iter_(niter-l), mean, covariance files and updated configuration files. 

To create joint marginal plots, simply run 

> plot_contours2d.pl -c /path/to/Sn/remap/config_pmc -n /path/to/Sn/remap 
/path/to/Lensing/remap 



New parameters, sampled from a flat or Gaussian distribution, can be added using add_par_f rom_prior . pi 



6.6. Analysis 

6.6.1. mvdens/mixjnvdens format utilities 



See Sect. A.3 for a description of the mvdens and mixjnvdens formats. 



fisher_tojneanvar .plj reads a mvdens file, inverts the covariance matrix and prints the mean 



and variance. 



corr_coeff . sh reads a mvdens or block matrix file and prints the correlation matrix of the 



covariance. 



diagjnvdens .pi replaces the covariance by its diagonal. 



add_par_to JTivdens . pi adds parameters to a mvdens file. Useful, if CosmoPMC is run with 
additional parameters, and the initial proposal is chosen from a previous run with the reduced 
parameter set. 
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6.6.2. PMC simulation/MCM chain utilities 



sampleZf ixpar . pi reads a sample file and fixes a parameter by cutting off" all points outside a 



given (narrow) range. 



6.6.3. PMC proposal diagnostics 

neff_proposal .plj calculates the effective number of components (eq.[5]l. It is the same quan- 
tity which is printed to the file enc. 



7. Using and modifying tlie code 

7.1. Modifying the existing code 

Note: Code to be used with MPI should not contain global variables and static variables. 

7.2. Creating a new module 

In this section, the steps required to add a new cosmology module to CosmoPMC are described. 

1. Create the directory newmodule and create (or copy) files with the necessary code to 
deal with the data and likelihood. Include files ("• .h) should be in newmodule/include, 
source files ( * . c) in newmodule/src. Edit the (or create a new) Makefile (in newmodule) 
and add the rules libnewmodule . so, libnewraodule . dylib and libnewmodule . a as 
well as the rule clean. 

2. In wrappers/include/types .h: 

Define a new data type by extending the enumeration data_t. Add the corresponding 
string (for identification of the module in the configuration file) in the macro sdata_t (i) , 
and increase Ndata_t by one. 

3. In wrappers/include/all wrappers. h: 
Add the line 

#include "newmodule .h" 

4. In tools/include/par .h: 

If necessary, add new parameter types (p jiewpararaeter) to enumeration par_t, add the 
corresponding identifier strings to the macro spar_t, and increase Npar_t by one. 

Optional: Add the parameter name and syntax for different programs (e.g. gnuplot, 
yorick, TgX) to bin/ spar . txt. 
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7. Using and modifying the code 



5. In wrappers/src/wrappers.c: 

Add the corresponding case to the 'switch' instruction in the function init_func_t. This 
function sets the data type. 

6. Create the files wrappers/include/newitiodule . h and wrappers/ src/newmodule . c. 
(Those files need to have different names than the files in newmodule/{src , include).) 
Write the following functions: 

a) init.functionjiewmodule 

b) read_from_config_newmodule 

c) init Jiewmodule 

d) likelijnewmodule (returning log L) 

e) special Jiewmodule (optional) 

f) print_newraodule (optional) 

To see what these functions are supposed to do, have a look at already existing modules, 
e.g. inbao.c. 

7. In Makefile .main: 

a) In the section "Additional directories", define the path to the new module's directory 

as 

NEWMODULE = $(COSMOPMC)/newmodule 

b) In the section "Libraries", define the library of the new module as 
LIBNEWMODULE = libnewraodule . $ (EXT) 

c) In the section "Combined cosmo include and linker flags", add the following flags: 
-1$ (NEWMODULE) /include to the variable IINCDIRS 

-L$ (NEWMODULE) to LLIBDIRS 
-Inewmodule to LLIBS. 

8. In exec/Makefile: 

Define the new rule: 

$(LIBNEWMODULE) : 

cd $ (NEWMODULE) && $(MAKE) $& 

(The second line has to start with a <TAB> and not with spaces.) 



9. Optional: Extend newdir_pmc . sh 
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7.3. Error passing system 

Most of the situations where an error occurs are intercepted by the program. In such a case, a 
variable '-err of type error* is set via the macros 

"err = addErrorCerror_type , "message", -err, LINE ); 

or 

"err = addErrorVA(error_type , "formatted message", '-err, 
__LINE__, VA_LIST); 

printing the current hne and function in the code, a message and the error type (negative integer). 
With 

testErrorRet(test , error_type, "message", '-err, LINE , 

return_value) ; 

or 

testErrorRetVA(test , error_type, "formatted message" , -'err, 
LINE , return_value, VA_LIST) ; 

a conditional error is produced if the (Boolean) expression test is true. The error is transported 
up the stack to the calling function with the macro 

forwardError(*err , LINE , return_value) ; 

Omit return.value in case of a void function. This can be used as diagnostics even for errors 
deep in the hierarchy of functions. 

During the calculation of the importance weights, any error is intercepted and the corresponding 
point does not contribute to the final sample. See Sect. |2] for more details. Therefore, in the 
routines which calculate the importance weights, the following is used: 

forwardErrorNoReturnC'-err , LINE , return_value) ; 

ParameterErrorVerbC'err , param, quiet, ndim) ; 

In case of an error, the first line forwards the error but does not return from the current rou- 
tine. The second line prints the ndim-dimensional parameter param to stderr (if quiet ! =1) and 
purges the error. 

To exit on an error, use 
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7. Using and modifying the code 



quitOnErrorC-err, LINE , FILE) 

This is usually done only from the main program. 

More macros and functions regarding error communication and handling can be found in the 
files errorlist .h, errorlist . c which are part of pmclib. 
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PMC references 



Introductory papers on PMC 



Cappe et al. 



Cappe et al. 



20041 



pOD8l 



Population Monte Carlo 

Adaptive importance sampling in general mixture classes 



Comparison of sampling methods including PMC 

Robert & Wraith (2009 1 Computational methods for Bayesian model choice 



Main papers on CosmoPMC 



Wraith et al. 



(20091 



Kilbinger et 



al. 



(20101 



Estimation of cosmological parameters using adaptive impor- 
tance sampling 

Bayesian model comparison in cosmology with Population 
Monte Carlo 



PMC applied to cosmological data 



Schrabback et al. 



Menard et al. 



(2010 



Benabed et al. 



Coupon et al. 



(120101) 



(20091 



(20121 



Kilbinger et al. 



Benjamin et al. 



Simpson et al. 



(20121 



(2012 



(20121 



Evidence of the accelerated expansion of the Universe from 

weak lensing tomography with COSMOS 

On the impact of intergalactic dust on cosmology with Type la 

supemovae 

TEASING: a fast and accurate approximation for the low mul- 
tipole likelihood of the cosmic microwave background tempera- 
ture 

Galaxy clustering in the CFHTLS-Wide: the changing relation- 
ship between galaxies and haloes since z. ~ 1.2-k 
CFHTLenS: Combined probe cosmological model comparison 
using 2D weak gravitational lensing 

CFHTLenS tomographic weak lensing: Quantifying accurate 
redshift distributions 

CFHTLenS: Testing the Laws of Gravity with Tomographic 
Weak Lensing and Redshift Space Distortions 



Other publications which use PMC 



Joachimi & Taylor 



Beaujean et al. 



(20I2I 



(|2011|) Forecasts of non-Gaussian parameter spaces using Box-Cox 
transformations 

Bayesian fit of exclusive b s€i decays: the standard model 
operator basis 
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A. File formats 



A.1. Data files 
A.1.1. Lensing 

For all lensdata.t types, the data format is the same. Each line contains the data for a given 
angular scale and (arbitrary many) redshift bin pair combinations. 

The angular scales are defined as follows. For lensformat = angle.center, the fist column 
contains the angular bin center in arc minutes. For the cases lensformat = angle_mean. 
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angle.wlinear and angle.wquadr, first two columns specify the lower and upper end of the 
angular bin. 

Following the angular- information are the data. For redshift bins, N^iN^ + l)/2 columns spec- 
ify all pair combinations {ij)i<j in lexical order, that is (1 1)(12)(13) . . . {\N^){22){23) . . . (N^Nz). 

Note that for lensdata = xipm the first No lines of the data file contain ^+ for Ng angular 
scales, the last Ne lines contain where the angular scales (first or first two columns) are 
identical in both halfs. 

The covariance matrix is in block format: It consists of N lines and N columns, where = 
NsN^iN^ + l)/2 is the length of the data. Usually, A^s is the number of measured angular scales, 
Ng, unless there is more than one data point per scale (e.g. for lensdata = xipm, = 2Ng). 

A matrix element equals {dtdj) - {di){dj), where di is the f-^ data point. In the counting 
over angular scale and redshift, the former varies faster than the latteij^ For example, with two 
redshift bins and three angular scales, the element C77 is the data variance for the redshift pair 
(11) and angular scale 6\ (starting counting at zero). Or, in other words, the covariance matrix 
consists of N^iNy + l)/2 block sub-matrices, each of size A^s x A^s- Each sub-matrix corresponds 
to one redshift bin combination. It is therefore easy to exclude some redshift bins, by (1) setting 
the diagonal of a sub-matrix to a very high value, and (2) setting the off-diagonal to zero (see 
the Nexclude parameter in the config file. Table 111. 



A.1.2. SNIa 

The SNIa data file in SN_SALT format starts with the following two lines: 

{aiNTRINSIC_DISPERSION double 
(aPECULIAR_VELOCITY double 

The peculiar velocity value is in units of km/s. This is followed by a list of supemovae, one 
object on each line as follows: 

2 2 7 

name z m s c <m><s><c'>< ms > < mc > < sc > 



A.1.3. BAO 

The BAO distance measures are modeled as Gaussian variables, the data files are in mvdens 
format (see Sect. A.3l. In the same file, following the mvdens data, there is a list of redshifts, 
corresponding to where the distances are measured. 



'This was wrongly stated here until version 1.01. 
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A.1.4. CMB 

The CMB data for WMAP are the ones released by the WMAP team. They are not included in 
CosMoPMC and can be obtained e.g. from the LAMBDA sit^ 

The SZ correction power spectrum file has two columns in each row containing £ and Cc, re- 
spectively. The first hne has to start with £ = 2. 

The CMB distance priors ( Komatsu et al.|2009 l are given in mvdens format. 



A.1.5. Redshift distribution 

The first line of a file describing a the redshift distribution for a redshift bin contains the type, 
see Sect. |4.7| 

# nofz 

This is followed by the list of parameter values, in the order given in Table [3] Each parameter 
value has to be in a new line, with the exeption of the histogram, nofz = single. In that case, 
the parameter lines are as follows: 

zo A^o 
Zl Ni 

Zn-l Nn-l 
Zn 

A^, is the number of galaxies in the bin [z,; z,+i]. The last hne denotes the upper limit of the last 
histogram bin z„ - Zmax> followed by a zero. For nofz = single, the file has to contain two 
identical lines with the value of zo in each line. 

A.2. Output file names 

The default names of all output files are defined in stdnames.h. Edit this file and to make 
clean; make to set user-defined file names. Note however that some of the pre-processing 
scripts expect the default names. 

A.3. Multi-variate Gaussian/Student-t (mvdens), mixture models (mixjnvdens) 

The mvdens file format is as follows. The first (header) line contains four integers: 



" http : / /lambda . gs£c . nasa . gov 
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p V B c. 

Here, p is the number of dimensions, v the degrees of freedom. For a multi-variate Gaussian, 
choose V - -\, and v > for Student- 1. B indicates the number of secondary diagonal of the 
covariance matrix which are updated during the PMC iterations. For most purposes, B can be 
set equal to p, which corresponds to the whole matrix being updated. Finally, c is 1 if the matrix 
is Cholesky-decomposed and otherwise. 

This is followed by p doubles indicating the mean, followed by p lines with p doubles each, 
giving the (symmetric) covariance matrix. 

Here is an example of a 5-dimensional multi-variate Gaussian (not Cholesky-decomposed): 
5-15® 

0.38559 -1.5238 19.338 1.3692 -2.4358 
0.0053677 -0.025608 0.00066748 -0.0011893 0.00087517 
-0.025608 0.16837 -0.0079163 0.0027364 -0.0035709 
0.00066748 -0.0079163 0.0011077 0.0010986 -0.00067815 
-0.0011893 0.0027364 0.0010986 0.016716 0.0026266 
0.00087517 -0.0035709 -0.00067815 0.0026266 0.014881 

The mixjtivdens format has two doubles as the header: 

D p 

where D is the number of components of the mixture and ndim the dimension. This is followed 
by D blocks specifying the weights Wd (doubles) and data md (in mvdens format) of the D 
multi-variate densities of the mixtures. 

m\ 

1712 
WD 

niD. 

The weights should be normalised, YJd=\ ^rf - 1- 

In many cases, an mvdens file indicates a parameter covariance matrix, for example to be used 
as Gaussian prior using the config file flag sprier. In some cases, the inverse covariance matrix 
is expected, as in the case of the Fisher matrix. 
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B. Syntax of all commands 



37 



All following scripts are located in $COSMOPMC/bin. All programs (executables) are located in 
$COSMOPMC/exec and linked from $COSMOPMC/bin after running make in SCOSMOPMC. 

• add deduced Jialomodel [TT|[22l 

Usage: add_deduced_halomodel [OPTIONS] PSIM [PAR_1 [PAR_2 [...]]] 
OPTIONS: 

-c CONFIG Configuration file (default: config_pmc) 

-o OUTNAME Ouput pmcsim name (default : psim+ded) 

PSIM pmc simulation file (pmcsim_iter) 

PAR_i String for deduced parameter #i . If not given, deduced 

parameters are read from the config file (default) 

• add_par_from_prior.pl 
Usage: add_par_from_prior.pl [OPTIONS] sample 

Adds a new random parameter to a PMC sample file, drawn under a distribution 
OPTIONS: 

Output sample file OUT (default: '<sample>.out' 
Prior distribution, DIST one of 'Flat' (default), 'Gauss' 
Prior arguments (white-spaced list if more than one) . For DIST 
Flat: ARG = 'min max' (defaut '-1 1') 
Gauss: ARG = 'mean sigma' 
Column COL of new parameter (default: last) 
Name string STR of new parameter 
This message 

add_par_to jnvdens . pi [37] 

add_par_to_mvdens.pl (MIX)MVDENS [OPTIONS] 

Adds a parameter to a (mix)mvdens file (e.g. Fisher matrix, PMC proposal) 
OPTIONS: 

Adds parameter in column and row COL (default : last column) 
Parameter mean VAL (default 0) 
Parameter variance VAL (default 
File is in 'mixmvdens' format 
File name 
This message 

bayes_factor.pl 



-o OUT 
-p DIST 
-P ARG 



-C COL 
-s STR 
-h 



-c COL 
-m VAL 
-V VAL 

-X 

FILE 
-h 



1) 
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Usage: bayes_factor.pl [OPTIONS] DIRl DIR2 

Calculates the Bayes factor between models. The corresponding 

evidence files (from PMC) have to be in the directories DIRl and DIR2 
OPTIONS: 

-i 'ITERl [ITER2] ' Use iteration ITERl for DIRl and ITER2 for DIR2 

(default: all iterations) 
-f 'EVIl [EVI2] ' Use files DIRl/EVIl and DIR2/EVI2 (default: 'evidence') 
-s Short output, last iteration only 

-1 Laplace approx. from Fisher matrix (denoted with iter=-l) 

-h This message 



cl_one_sided [10] 
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Usage: cl_one_sided [OPTIONS] sample 
OPTIONS: 



-c CONFIG 
-i INDEX 
-d DIR 
-V VALUE 
-w WHICH 

sample 

The options 



Configuration file (default: config_pmc) 
Parameter index 
Direction (DIR=+1,-1) 
Starting value 

WHICH=0: 68%, 95%, 99. 7% c.l. (default) 
WHICH=1: 68%, 90%, 95% c.l. 
PMC sample file 
INDEX, -d DIR and -v VALUE are required 



config_pmc_tojnax_and_fish.pl |7 26 



Usage: config_pmc_to_max_and_fish.pl [OPTIONS] 
OPTIONS: 

Create config file for maximum search (max_post) 
Create config file for Fisher matrix (go_fishing) 
Input PMC config file CONFIG (default: 'config_pmc') 
Random starting point (for maximum search) 
Fiducial starting point FID. FID is a white-space 

separated list in quotes, e.g. '0.25 0.75' 
Fidcucial parameter from FILE (e.g. 'maxlogP') 
Tolerance for maximum- search (default: 0.01) 
Calculate only diagonal of Fisher matrix (go_fishing) 
This message 
'-F' is obligatory 



-M 
-F 

-c CONFIG 
-r 

-f FID 

-p FILE 
-t TOLERANCE 
-d 
-h 

One of '-M' or 

The default starting point for maximum search is (max-min)/2 
For Fisher matrix ('-F'), a fiducial parameter has to be indicated with 
or '-p FILE' 

corr_coeff . sh [37] 

Usage: corr_coeff filename [mvdens | block] 
cosmojncmc 

Usage: cosmo_mcmc [OPTIONS] 
OPTIONS : 



-f FID' 



-c CONFIG 
-s SEED 

-h 

cosmo_pmc 

Usage: cosmo_pmc [OPTIONS] 
OPTIONS : 



Configuration file (default: config_mcmc) 
Use SEED for random number generator. If SEED=-1 (default) 

the current time is used as seed. 
This message 
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-c CONFIG 
-s SEED 

-q 
-h 



Configuration file (default: 'config_pmc') 

Use SEED for random number generator. If SEED=-1 (default) 

the current time is used as seed. 
Quiet mode 
This message 



cosmo_pmc.pl Blu 35 



Usage: cosmo_pmc.pl [OPTIONS] 
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OPTIONS : 
-n NCPU 
-c CONFIG 
-f FID 



-m [c|a] 
-d 



SEED 
[M|F] 

[y|n] 

PATH 



-e 

-p PRO 
-M MULT 
-0 OPT 

-q 

-h 



Run PMC in parallel on NPCU cpus using 'mpirun' (default: 1) 
Configuration file for PMC (default: config_pinc) 
Fiducial starting point FID. FID is a white-space 

separated list in quotes, e.g. '0.25 0.75' 
Random starting point for maximum search 

(default: (max-min)/2) 
Maximum- search method: 'c' (eg), 'a' (amoeba) 
Calculate only diagonal of Fisher matrix 

Do not force Fisher matrix F to be positiv. If F is negative, 

script exits with an error 
Adaptive numerical differentiation for Fisher matrix 
Use SEED for random number generator. If SEED=-1 (default) 

the current time is used as seed. 
Stops after maximum search ('M') or Fisher matrix ('F') 
Default answer to all questions on stdin 
Use PATH as CosmoPMC directory (default : environment 

variable SCOSMOPMC) 
Create 'essential' plots 

Plotting scripts: 'y' (yorick; default), 'R' (R) or 'n' (none) 
Combinations of letters are possible, e.g. 'yR' 

Output sample MULT times input (default 1) . 
Valid if plotting script is 'R' 

Pass options OPT to 'plot_contour2d.pl' 

Quiet mode 

This message 



diag_mvdens.pl [37] 

Usage: diag_mvdens.pl IN 

Prints the mvdens file 'IN' with the covariance replaced by its diagonal. 

essential_cosmo_pmc_run . pi 
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-c 

-p 
-k 

-V 

-h 



CONFIG 
PATH 



Usage: essential_cosmo_pmc.pl [OPTIONS] 
OPTIONS: 

Uses config file CONFIG (default: 'config_pmc') 
Use PATH as CosmoPMC directory (default : environment 

variable SCOSMOPHC) 
Keep temporary files 
Verbose 
This message 



evidence.pl [36] 



Usage: evidence.pl [OPTIONS] SAMPLE 
OPTIONS: 

-h This message 

SAMPLE PMC sample file 

evidence.list . pi 

Usage: evidence_list.pl [OPTIONS] DIRl [DIR2 [. 
OPTIONS : 

-r N Subtract log(E) from DIRN (default: no subtraction) 



36 



..]] 
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For N=-l subtract log(E_min) 
-k KEY Use KEY (string list) instead of 

directory names (default) 
-s SEP Use SEP as input separator for KEY list 

-S SEP Use SEP as output separator 

(default for both: white-space) 
-n Write number of model parameters 

-L Use Laplace approximation (reading file ' evidence_fisher ' ) 

-h This message 

fisher_to_meanvar.pl |37] 

fisher_to_meanvar.pl [OPTIONS] file 
OPTIONS: 

-n No inverse 

-m Marginal errors (don't invert matrix) 

-X mixmvdens format (default: mvdens format) 

-k Keep temporary file 'fishtmp.i' 

-h This message 

Options '-m' and '-n' exclude each other 

get_spar.pl 

Usage: get_spar.pl [OPTIONS] LANG [PARI [PAR2 [...]]] 
OPTIONS : 

-c CONFIG Configuration file ONFIG (default 'config_pmc') 

-i INDEX Returns only par [INDEX] 

-P PATH Use PATH as CosmoPMC directory (default: environment 

variable SCOSHOPMC) 

-p Print 'p<i> for unknown parameters instead of input string 

LANG One of 'yorick' , 'gnuplot', 'TeX', 'R' . 

More languages can be defined in spar.txt 
PARI . . . Prameter strings 

go_fishing |7j[l2j[T2] 

Usage: go_fishing [OPTIONS] 
OPTIONS : 

-c CONFIG Configuration file (default: config_fish) 

-a Adaptive numerical differentiation (default: fixed difference) 

-f Force positive Fisher matrix 

-q Quiet mode 

-h This message 
Run in parallel on NP cpu's: 'mpirun -np NP go_fishing [OPTIONS] 



haloplot 22 



Usage: haloplot loglO(M_min) log 18 (Ml) logl8(M8) sigma_log_M alpha_halo 

halomodel.par [OPTIONS] 
Outputs HOD-derived quantities 
OPTIONS: 

-o OUT Output file name 

-t TYPE Output type, TYPE in [wtheta, wp, xi , xihalo, deltaSig, 

nofm, halo, pk] , default: wtheta 
-nbins Number of bins 
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-range Range (linear scale): min,max 

-z Z Used fixed redshift Z (no w(theta) output) 

-Mhalo loglOM log 10 (Halo mass) for deltaSig and xihalo (in M_sol/h) 

-c CONFIG PMC config file, to calculate chi"2 

-h This message 

M_min, Ml and MO are in units of M_{sol}/h. 

histograms_sample [TT][TT] 

Usage: histograms_sample [OPTIONS] sample 
OPTIONS: 

-c CONFIG Configuration file (default: 

-1 Only Id histograms 

-2 Only 2d histograms 

sample PMC sample file 

-h This message 

importance_sample 



config_pmc) 
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-c CONFIG 
-o OUTSAMPLE 

-q 

-h 

INSAMPLE 



-c CONFIG 
-m [c|a|n] 

-t 

-s SEED 



config_pmc) 
'insample.out') 



Usage: importance_sample [OPTIONS] INSAMPLE 
Performs an importance run on a PMC sample . Run in 

parallel with MPI (use mpirun) 
OPTIONS: 

Configuration file (default: 
Output sample name (default: 
Quiet mode 
This message 
Input sample name 

max.post |7j[T2|[T2] 

Usage: max_post [OPTIONS] 
OPTIONS: 

Configuration file (default: config_max) 
Maximum-search method: 'a' (amoeba, default), 'c' (eg), 

'n' (none; print posterior for fiducial parameter and exit) 
Test maximum at the end 

Use SEED for random number generator. If SEED=-1 (default) 

the current time is used as seed. 
Prints the maximum-posterior model to the file 'modeljiaxlog' 
Quiet mode 
This message 



-P 

-q 

-h 



mean2eps.pl 35 



Usage: mean2eps.pl [OPTIONS] MEAN 
OPTIONS: 

MEAN File containing mean and confidence levels (output of 

'cosmo_pmc' or 'histograiiis_sample' 
-c CONFIG Uses config file CONFIG (default: 'config_pmc') 

-P PATH Use PATH as CosmoPMC directory (default : environment 

variable SCOSMOPMC) 
-o BASE Outname BASE (default: <MEAN>) 

-V Verbose 
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-h 



This message 



meanvar2tab.pl 35 



Usage: meanvar2tab.pl [OPTIONS] file [file2 [...]] 



Options: 
-s {123} 
PREC 



-P 
-e 
-c 
-t 



CONFIG 
TITLE 



-S SEP 
-P PATH 



68% (1), 95% (2) or 99.7% (3) errors (default = 1) 

Output with PREC digits ('%PREC' format string) 

Error(s) written to PREC significant digits (use -p PREC) 

Uses config file CONFIG (default: 'config_pmc') 

Title (table heading) TITLE is string list with entries according 

to the number of input files 

Use SEP as input separator for TITLE list (default: white space) 
Use PATH as CosmoPMC directory (default: environment 

variable SCOSMOPMC) 
This message 



meanvar_sample 10 



Usage: meanvar_sample [OPTIONS] sample 
OPTIONS : 



-c CONFIG 
-w 
-C 
-E 

-h 

sample 
nef f _proposal . pi 
Usage: neff_proposal.pl PROP 



Configuration file (default: config_pmc) 

Ignore weights (default: weights=first column of sample file) 
Write covariance and inverse covariance to files 
Output evidence 
This message 
PMC sample file 

T0|[38 



Calculates the effective number of components for the mix_mvdens file 'PROP' 
newdir_pmc . sh [5|[39 



Usage : newdir_pmc . sh [DIR] 

Directory DIR (default: read on input) is created. 

Links are set to data files in \$COSMOPMC/data. 

Parameter files are copied on request from \$COSMOPMC/par_files . 

plot confidence. R [6| [TT|[35 



Usage: plot_confidence.R [options] 



Options: 
-h, --help 

Show this help message and exit 
-N NGRID, — Ngrid=NGRID 

Number of grid points for smoothing (kde2d) (default 180) . Use <=30 for 
fast-but-dirty plots 

-g GSMOOTH, --gsmooth=GSMOOTH 

Smoothing kernel width, with respect to box size (default 38). In case of more 



54 



B. Syntax of all commands 



than one sample, use list separated with '_' for more than value 
-S, --solid 

All contours with solid lines 

-w WIDTH, --width=WIDTH 
Line width (default 1) 

-k, --with_keys 
Add key to plots 

-K KEYSTRING, --keystring=KEYSTRING 
Key strings (separate items with '_') 

-L, --no_key_line 

Do not add a line to the keys in the legend 

-c CONFIG, --config=CONFIG 

Config file (default 'config_pmc') 

-t TITLE, --title=TITLE 

Title string for each panel (default empty) 

-i INDEX_I, --index_i=INDEX_I 

Only create plots with i-th parameter on x-axis 

-j INDEX_J, --index_j=INDEX_J 

Only create plots with j-th parameter on y-axis 

-s SIGMA, --sigma=SIGMA 

Plot SIGMA confidence levels (default 3) 

-F COLOR_SCHEME , --color_scheme=COLOR_SCHEME 
Color scheme (0, 1; default 8) 



plot contour2d.pl [6| [TT] [TT] [35 

Usage: plot_contour2d.pl [OPTIONS] [DIRl [DIR2 [...]]] 



OPTIONS : 

-i NITER Number of iterations (needed if do_proposal=2) 

-c CONFIG_FILE Configuration file (default: in order config_mcmc, config_pmc) 

-t TITLE Title string for each panel (default empty) 

-T TITLE Title string for all_contour2d. {eps |pdf} (default empty) 

-n No shade 

-w WIDTH Line width WIDTH (default 4) 

-1 OPT Add Id posterior plots. OPT can contain the following letters: 

m Plot line at mean position 
123 Plot line at 68%, 95%, 99. 7 density 
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-s 

-s N 

-r 

-g FACTOR 



-G FACTOR 
-C 

-N NORM 



-F NUM 
-k 

-K "KEYl 
-y FS 
-o FORMAT 
-b 

-m PAR 
-P PATH 

-q 
-h 

DIRl . . . 



[KEY2 



t Write mean and 68% confidence intervals as text 
(use with 'm' and ' 1 ' 

n None of the above 
All contours with solid lines 
Outermost level is N sigma 

Aspect ratio=l, changes plot limits such that dx=dy 
Gaussian smoothing of 2d-histograms with variance 

box-width/ 1 FACTOR I . If FACTOR is negativ, plots 

unsmoothed histogram in addition (use with '-n'). 

Note: For multiple contours, use a list of values "gl g2 ..." 
Gaussian smoothing of Id-histograms (default: 2d factor) 
Use covariance (file covar.fin) for Gaussian smoothing 
Normalisation of Id posterior 

'm' Maximum = 1 (default) 

'i' Integral over posterior = 1 
Color scheme, NUM=8,1,2 
Add key to plots 

[...]]" Key strings (default: directory names) 
Font size FS (default 24) 

Output file format, FORMAT=eps |pdf (default: eps) 
Writes the chi2 files in block format 

Plots a mark at position PAR (e.g. best-fit). PAR is white-space 
separated list (use quotes or '\ ', e.g. '8.3 0.8') 

Use PATH as CosmoPMC root directory (default: environment 
variable SCOSMOPMC) 

Run quietly, no verbose 

This message 

List of directories containing histogram files (chi2_*_*) 
Default: DIRl = ' . ' 



proposal_mean.pl [T0l[36 



Usage: proposal_mean.pl [OPTIONS] 
OPTIONS: 



-d DIR 



CONFIG 



PATH 



Directory DIR containing the sub-directories 'iter_*' 

with the proposal files (default '.') 
Configuration file CONFIG (default 'DIR/config_pmc') 
No plotting, only creates '.gnu' file 
X- and y-axes inverted 
X- and y-labels on top/right 

Use PATH as CosmoPMC root directory (default: environment 

variable SCOSMOPMC) 
This message 



proposal_var.pl [10} [36 



Usage: proposal_var.pl [OPTIONS] 
OPTIONS: 

-d DIR Directory DIR containing the sub-directories 'iter_*' 

with the proposal files (default '.') 
-c CONFIG Configuration file CONFIG (default 'DIR/config_pmc') 

-P PATH Use PATH as CosmoPMC root directory (default: environment 
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36 



-h 

remap . sh 

Usage: remap. sh [OPTIONS] 
OPTIONS: 



variable SCOSMOPMC) 
This message 



-c CONFIG 
-i INPUT 
-s PMCSIM 
-o OUTPUT 
-r REMAP 
-n NPAR 
-d N_DED 
-h 



Input PMC configuration file (default ' ./config_pmc') 

Input directory INPUT (default '.') 

Sample/PMC simulation file PMCSIM 

Output directory OUTPUT (default ' . /remap ' ) 

Remap file REMAP (default '. /remap . dat ' ) 

Number of parameters NPAR (default: read from remap file) 

Number of deduced parameters N_DED (default: 8) 

This message 



sample2fixpar.pl [38] 

Usage: sample2fixpar.pl SAMPLE_IN COL MIN MAX 



SAMPLE_IN 
COL 



MIN, MAX 
tabZtex . pi 
Usage: tab2tex.pl [OPTIONS] file 



Input sample (PMC simulation or MCM chain) 
Column number of fixed parameter 

(Note that par #i is in column i+2) 
Minimum and maximum values for fixed parameter 
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OPTIONS: 

-a 
-b 

-s STRETCH 
-m 

-1 MODE 



-L MODE 



-h 



Produce tex array, not tex table 

Bare output , no table/array header 

Set arraystretch to STRETCH 

Add '$' around entries (tex inline math mode) 

Print vertical lines between rows according to MODE; 

a all lines (default) 

n no lines 

h header lines 
Print horizontal lines between columns according to Mode: 

a all lines (default) 

n no lines 
This message 



test_suite_cosmo_pmc . pi 

Usage: test_suite_cosmo_pmc.pl [OPTIONS] 
OPTIONS: 

-r Do PMC test runs 

-R Only do PMC test runs 

-n NCPU Run PMC in parallel on NCPU cpus (default: 1) 

-c Include CMB tests 

-P PATH Use PATH as CosmoPMC root directory (default: environment 

variable SCOSMOPMC) 

-s Short, without time-taking PMC runs (e.g. Lensing/COSMOS-SlS) 

-k Keep temporary files 

-X Clean previous run and exit 
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start 




read chain. fin 



f 



no 



alculate marginals, 
create histograms 




remove bum-in, 
decorrelate 



no 



write chain. fin 




read chain.pre 



t-unMCMC,i^ 
rite chain.acc ] 



initial proposal 




Figure 2: Flow chart of the MCMC implementation. 



-V 

-h 



Verbose 
This message 



C. MCMC 



We provide a Metropolis-Hastings Monte-Carlo Markov Chain sampler, which is included in 



the CosMoPMC package. This MCMC implementation has been used in Wraith et al. (20091 in 
the comparison with PMC. In the following, we briefly describe our MCMC program. 
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Table 13: MCMC section of the configuration file 



nchain 
ncov 
fburnin 
ndecorr 

fudge 
sinitial 



boxdiv|^ 
sstart 



integer 
integer 
double 
double 

double 
string 



double 
string 



Chain Length 

Interval between updates of the proposal covariance 
Bum-in phase are the first ncovxncor points 
De-correlation (thinning-out): one in ndec points is kept 
in the final chain 

Proposal covariance is multiplied by fudge^/n_par 
Initial proposal type, one of Fisher_inv, Fisher, 
Fisher, previous, Hessian, Hessian_diag , 
diag. 

Diagonal of proposal covariance is (max-min)/boxdiv 
Starting point type, one of ran , fid, min, max, nul 



npar doubles Starting parameter 



Histogram section 



nbinhist integer 



Number of density histogram bins 



"only if sinitial = diag 
*only if sstart = fid 



C.1. MCMC configuration file 
C.2. Proposal and starting point 

The proposal for the Metropolis-Hastings algorithm is a multi-variate Gaussian distribution. 
After choosing an initial proposal, a new proposal can optionally be re-calculated after a number 
of ncov (accepted) steps. The covariance of this new proposal is the chain covariance from steps 
up to this point. This proposal is then updated after each ncov accepted steps using all previous 
accepted points. 

There are several options for the initial proposal: 

1. sinitial = diag A diagonal covariance with width being a fraction of the box size. 

2. sinitial = FisherTheHessianatagivenpointinparameter space. If this point is the 
maximum-likelihood point, the Hessian corresponds to the Fisher matrix. 

3. sinitial = Fisher_inv The inverse Hessian/Fisher matrix, e.g. the covariance from a 
previous chain. This can be useful for ill-conditioned matrices which are difficult to invert 
numerically. 

4. sinitial = previous A proposal read from a file, e.g. from a previous MCMC run. 

The starting point is either chosen randomly or specified in the config file. The second case 
might be convenient if the prior volume is very large and a very long burn-in phase is to be 
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avoided. For example, the ML-point or best-fit value from a previous experiment can be chosen 



Dunkley et aLl ( l2009| ). 



C.3. Output files 



The MCMC output files have the same format as their PMC counterparts (see Sect. 3.3.2). 

A complete run of cosmo_mcmc produces three files containing the points of the Markov chain: 

1. chain. all containing all, accepted and rejected, sample points. This is the only chain 
file will not be read or used in subsequent calls of cosmo jncmc. 

2. chain . acc containing the accepted points. 

3. chain . fin containing the accepted points after removal of the bum-in phase and after de- 
correlating (thinning-out) the chain. The results produced by cosmojticmc (mean, errors, 
histograms, covariance) are based on this file. 

The chains are ASCII-files, in the same format as the PMC sample files. All weights are 1, and 
the second column contains the log-likelihood (only in chain . all. 

The parameter mean and confidence intervals are printed to the file mean. The names of files 
containing the histograms and parameter covariances are the same as for PMC. 

C.4. Diagnostics 

In general it is not straight-forward to diagnose an MCM chain. There exists tests but no formal 
proofs for convergence (e.g. Gellman-Rubin), which in addition require very long or multiple 
chains. We have not implemented such tests in the code. However, there are a few (rather 
hand-waving) diagnostic tools to check the reliability of an MCMC run. 

Firstly, the acceptance rate 77 should be in the range between 15% and 25%. A larger 77 most 
probably corresponds to a chain which stayed mainly in the high-density region and strongly 
under-sampled the lower-density posterior regions. In that case the error bars will be underesti- 
mated. A very small 77 means probably an under-sampling of the posterior since only few points 
are accepted. However, this need not cause a bias for the parameters and errors if the chain has 
been run long enough. 

C.5. Resuming an interrupted run 

Sometimes a MCMC run is interrupted before finishing, or one wishes a previous run to be 
extended, for example because its convergence is doubted. The MCMC program allows to read 
in and extend a previous chain. To that end, rename the file chain . acc into chain . pre. The 
proposal for the resumed run can but need not be calculated from the previous chain (to be 
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controlled in the config file, see Sect C.2i. In the config file, the number of desired sample 



points has to be larger than the previous chain. 
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